Loading Libraries

library(tidyverse)
library(sf)
library(leaflet)
library(htmltools)
library(htmlwidgets)
library(raster)
library(gstat)
library(spatial)

Loading Neighborhood Data

Data from Analyze Boston and the BPDA 2019 Annual Neighborhood Profiles (http://www.bostonplans.org/getattachment/f719d8d1-9422-4ffa-8d11-d042dd3eb37b)

boston_nhoods <- st_read("https://opendata.arcgis.com/datasets/3525b0ee6e6b427f9aab5d0a1d0a1a28_0.geojson", quiet = TRUE) %>%
  dplyr::select(Name)

leaflet(boston_nhoods) %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addPolygons(highlightOptions = highlightOptions(fillColor = "yellow", 
                                                  fillOpacity = 1),
              label = ~Name, 
              weight = 1) 

Cloropleth Map

boston_nhoods <- boston_nhoods %>%
  mutate(med_income = c(68209, 76968, 32103, 38235, 65090, 65090 ,65090, 82965, 25937, 77161, 93033, 51549, 91067, 90694, 93033, 65090, 34483, 50110, 79424, 65260, 43256, 47200, 111518, 77223, 88469, NA))

boston_nhoods2 <- drop_na(boston_nhoods)

boston_nhoods2$label <- paste("<b>",
                               "<p style='color:Black;'>", 
                               boston_nhoods2$Name,"</b>","</p>",
                               "<b>", "Median Income: ",
                               ifelse(boston_nhoods2$med_income > 62021,
                                      "<span style = 'color:Green;'>", 
                                      "<span style = 'color:Red;'>"), "",
                               "$",
                                prettyNum(boston_nhoods2$med_income, 
                                             digits = 6, big.mark = ","),"</b>", "</span>") %>%
  lapply(htmltools::HTML)

bins <- seq(min(boston_nhoods2$med_income),
            max(boston_nhoods2$med_income))
pal <- colorNumeric("viridis", 
                    domain = boston_nhoods2$med_income,
                    na.color = "#00000000")

leaflet(boston_nhoods2) %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addPolygons(highlightOptions = highlightOptions(fillOpacity = 1),
              label = ~label,
              fillColor = ~pal(med_income),
              weight = 1, color = "black") %>% 
  
addControl("Median Neighborhood Income, Boston MA", position = "topright") %>%

  
addLegend(pal = pal, 
            values = ~med_income,
            bins = 4,
            opacity = 0.7, title = "Median Income (2017)",
            position = "topright")%>%

addLegend(position = "bottomright",
            colors = c("green","red"),
            labels = c("Neighborhood median income is higher than Boston's",
                       "Neighborhood median income is lower than Boston's"),
            title = "Census Tract Pop-up Text")  %>%
  
    addScaleBar(position = "bottomleft") %>%


setView(lng = -71.037406, lat = 42.339542, zoom = 11) %>%

setMaxBounds( lng1 = -71,
                lat1 = 42,
                lng2 = -71.1,
                lat2 = 42.5)

Area Centroid Map

MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs "

WGS84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

points <- st_centroid(
  st_transform(boston_nhoods, crs = MA_state_plane)) %>%
  st_transform(WGS84)

map2 <- leaflet(points) %>%
  addControl("Median Neighborhood Income, Boston MA", position = "topright") %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addCircles(highlightOptions = highlightOptions(fillOpacity = 1),
             fillColor = ~pal(med_income),
             stroke = FALSE, 
             radius = 100, 
             fillOpacity = 1) %>%
  
  addLegend(pal = pal, 
            values = ~med_income,
            bins = 4,
            opacity = 0.7, title = "Median Income (2017)",
            position = "topright")  %>%
  
  addScaleBar(position = "bottomleft") %>%


  setView(lng = -71.037406, lat = 42.339542, zoom = 11.48) %>%

  setMaxBounds( lng1 = -71,
                lat1 = 42,
                lng2 = -71.1,
                lat2 = 42.5)

map2
saveWidget(map2, file = "Boston_Income_Map2.html")

Continuous Surface Map

nhoods_inter <- points %>%
  st_transform(MA_state_plane) %>%
  as_Spatial()

nhoods_inter <- boston_nhoods2 %>%
  st_transform(MA_state_plane) %>%
  as_Spatial()

boston_raster <- raster(nhoods_inter, res=10)

gs <- gstat(formula=med_income~1, locations=nhoods_inter)
idw_interp <- interpolate(boston_raster, gs)
## [inverse distance weighted interpolation]
idw_interp_clip <- mask(idw_interp, nhoods_inter)
map3 <- leaflet(points) %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addRasterImage(idw_interp_clip, colors = pal, opacity = 0.8) %>% 
  addControl("Median Neighborhood Income, Boston MA", position = "topright") %>%
  addLegend(pal = pal, 
            values = ~med_income,
            bins = 4,
            opacity = 0.7, title = "Median Income (2017)",
            position = "topright") %>%

  addScaleBar(position = "bottomleft") %>%


  setView(lng = -71.037406, lat = 42.339542, zoom = 11.48) %>%

  setMaxBounds( lng1 = -71,
                lat1 = 42,
                lng2 = -71.1,
                lat2 = 42.5)

map3
saveWidget(map3, file = "Boston_Income_Map3.html")
boston_nhoods <- boston_nhoods %>%
  mutate(med_income = c(68209, 76968, 32103, 38235, 65090, 65090 ,65090, 82965, 25937, 77161, 93033, 51549, 91067, 90694, 93033, 65090, 34483, 50110, 79424, 65260, 43256, 47200, 111518, 77223, 88469, NA))

boston_nhoods2 <- drop_na(boston_nhoods)

boston_nhoods2$label <- 
   paste("<b>",  
         "<p style='color:Black;'>",
          boston_nhoods2$Name, "</b>","</p>",
          "Median Income:",
          ifelse( boston_nhoods2$med_income==0,
                                      "<span style = 'color:Gray;'> Not Available",
                                      ifelse(boston_nhoods2$med_income >
                                               62021,
                                             "<span style = 'color:Green;'>",
                                             "<span style = 'color:Red;'>")),
                               ifelse(boston_nhoods2$med_income == 0, "",
                               "$"), ifelse(boston_nhoods2$med_income == 0, "",
                                 prettyNum( boston_nhoods2$med_income, 
                                             digits = 6, big.mark = ",")),
                               "</span></b></p>") %>%
  
   lapply(htmltools::HTML)


bins <- seq(min(boston_nhoods2$med_income),
            max(boston_nhoods2$med_income))
pal <- colorNumeric("viridis", 
                    domain = boston_nhoods2$med_income,
                    na.color = "#00000000")

map4 <- leaflet(boston_nhoods2) %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addControl("Median Neighborhood Income, Boston MA", position = "topright") %>%
  addPolygons(highlightOptions = highlightOptions(fillOpacity = 1),
              label = ~label,
              fillColor = ~pal(med_income),
              weight = 1, color = "black") %>% 
  addLegend(pal = pal, 
            values = ~med_income,
            bins = 4,
            opacity = 0.7, title = "Median Income (2017)",
            position = "topright") %>%
  addLegend(position = "bottomright",
            colors = c("green","red"),
            labels = c("Higher than Boston Median Income",
                       "Lower than Boston Median Income"),
            title = "Census Tract Pop-up Text") %>%  

 addScaleBar(position = "bottomleft") %>%


  setView(lng = -71.037406, lat = 42.339542, zoom = 11) %>%

  setMaxBounds( lng1 = -71,
                lat1 = 42,
                lng2 = -71.1,
                lat2 = 42.5)

map4
saveWidget(map4, file = "Boston_Income_Map4.html")

Map Comparisons

Most Informative Map

In my opinion, the most informative map is the choloropleth map because it allows users to scroll over the neighborhoods and clearly understand the median income through the pop that lists it versus requiring users to to look at the color shade and estimate the median income for the location by comparing that shade to the key.

Most Interesting Map

In my opinion, the most interesting map is the continuous map because it is visually the most interesting to look at and it adds new information that the other maps do not through interpolation. This map also highlights how neighborhood boundaries, while sometimes driven by spatial factors are often driven by politics. While in the chloropleth map we see stark comparisons between neighborhoods, in the interpolated map we see how median income could vary between two neighborhood boundaries and move beyond political lines.

Most appropriate to the data

I think the most appropriate for the data is the chloropleth map because it most accurately represents the data collected. The continuous map and the area centroid map assume that the median income value represents the income at the center of the neighborhood and ripples out continuously, however, without more fine-grained data, I think visualizing the data this way could be inaccurate and misleading.

Best Map

For the reasons listed in the “most appropriate” section, I think the choloropleth map is the best map out of the three. I think it is more legible than the centroid map and more accurate than the interpolation map, making it, in my opinion, the all around best map.